home *** CD-ROM | disk | FTP | other *** search
/ Aminet 30 / Aminet 30 (1999)(Schatztruhe)[!][Apr 1999].iso / Aminet / gfx / misc / gnuplot-3.7src.lha / gnuplot-3.7src / gnuplot-3.7.lha / gnuplot-3.7 / plot3d.c < prev    next >
C/C++ Source or Header  |  1998-12-10  |  48KB  |  1,577 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: plot3d.c,v 1.36 1998/06/18 14:55:14 ddenholm Exp $";
  3. #endif
  4.  
  5. /* GNUPLOT - plot3d.c */
  6.  
  7. /*[
  8.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  9.  *
  10.  * Permission to use, copy, and distribute this software and its
  11.  * documentation for any purpose with or without fee is hereby granted,
  12.  * provided that the above copyright notice appear in all copies and
  13.  * that both that copyright notice and this permission notice appear
  14.  * in supporting documentation.
  15.  *
  16.  * Permission to modify the software is granted, but not the right to
  17.  * distribute the complete modified source code.  Modifications are to
  18.  * be distributed as patches to the released version.  Permission to
  19.  * distribute binaries produced by compiling modified sources is granted,
  20.  * provided you
  21.  *   1. distribute the corresponding source modifications from the
  22.  *    released version in the form of a patch file along with the binaries,
  23.  *   2. add special version identification to distinguish your version
  24.  *    in addition to the base release version number,
  25.  *   3. provide your name and address as the primary contact for the
  26.  *    support of your modified version, and
  27.  *   4. retain our contact information in regard to use of the base
  28.  *    software.
  29.  * Permission to distribute the released version of the source code along
  30.  * with corresponding source modifications in the form of a patch file is
  31.  * granted with same provisions 2 through 4 for binary distributions.
  32.  *
  33.  * This software is provided "as is" without express or implied warranty
  34.  * to the extent permitted by applicable law.
  35. ]*/
  36.  
  37. #include "plot.h"
  38. #include "setshow.h"
  39. #include "binary.h"
  40.  
  41. #ifndef _Windows
  42. # include "help.h"
  43. #endif
  44.  
  45. #ifndef STDOUT
  46. #define STDOUT 1
  47. #endif
  48.  
  49. /* static prototypes */
  50.  
  51. static void get_3ddata __PROTO((struct surface_points * this_plot));
  52. static void print_3dtable __PROTO((int pcount));
  53. static void eval_3dplots __PROTO((void));
  54. static void grid_nongrid_data __PROTO((struct surface_points * this_plot));
  55. static void parametric_3dfixup __PROTO((struct surface_points * start_plot, int *plot_num));
  56.  
  57. /* the curves/surfaces of the plot */
  58. struct surface_points *first_3dplot = NULL;
  59. static struct udft_entry plot_func;
  60.  
  61.  
  62. extern struct udft_entry *dummy_func;
  63. extern int datatype[];
  64. extern char timefmt[];
  65. extern TBOOLEAN is_3d_plot;
  66. extern int plot_token;
  67.  
  68. /* in order to support multiple axes, and to
  69.  * simplify ranging in parametric plots, we use
  70.  * arrays to store some things. For 2d plots,
  71.  * elements are  y1 = 0 x1 = 1 y2 = 2 x2 = 3
  72.  * for 3d,  z = 0, x = 1, y = 2
  73.  * these are given symbolic names in plot.h
  74.  */
  75.  
  76. extern double min_array[AXIS_ARRAY_SIZE], max_array[AXIS_ARRAY_SIZE];
  77. extern int auto_array[AXIS_ARRAY_SIZE];
  78. extern TBOOLEAN log_array[AXIS_ARRAY_SIZE];
  79. extern double base_array[AXIS_ARRAY_SIZE];
  80. extern double log_base_array[AXIS_ARRAY_SIZE];
  81.  
  82. /* some file-wide variables to store which axis we are using */
  83. static int x_axis, y_axis, z_axis;
  84.  
  85.  
  86. /* Deleted from setshow.h and renamed */
  87. extern FILE *gpoutfile;
  88.  
  89. /* info from datafile module */
  90. extern int df_datum;
  91. extern int df_line_number;
  92. extern int df_no_use_specs;
  93. extern int df_eof;
  94. extern int df_timecol[];
  95. extern TBOOLEAN df_matrix;
  96.  
  97. #define Inc_c_token if (++c_token >= num_tokens)    \
  98.                         int_error ("Syntax error", c_token);
  99.  
  100. /* From plot2d.c */
  101. extern int reverse_range[];
  102.  
  103. /*
  104.  * IMHO, code is getting too cluttered with repeated chunks of
  105.  * code. Some macros to simplify, I hope.
  106.  *
  107.  * do { } while(0) is comp.lang.c recommendation for complex macros
  108.  * also means that break can be specified as an action, and it will
  109.  * 
  110.  */
  111.  
  112. /*  copy scalar data to arrays
  113.  * optimiser should optimise infinite away
  114.  * dont know we have to support ranges [10:-10] - lets reverse
  115.  * it for now, then fix it at the end.
  116.  */
  117. #define INIT_ARRAYS(axis, min, max, auto, is_log, base, log_base, infinite) \
  118. do{if ((auto_array[axis] = auto) == 0 && max<min) {\
  119.     min_array[axis] = max;\
  120.    max_array[axis] = min; /* we will fix later */ \
  121.  } else { \
  122.     min_array[axis] = (infinite && (auto&1)) ? VERYLARGE : min; \
  123.     max_array[axis] = (infinite && (auto&2)) ? -VERYLARGE : max; \
  124.  } \
  125.  log_array[axis] = is_log; base_array[axis] = base; log_base_array[axis] = log_base;\
  126. }while(0)
  127.  
  128. /* handle reversed ranges */
  129. #define CHECK_REVERSE(axis) \
  130. do{\
  131.  if (auto_array[axis] == 0 && max_array[axis] < min_array[axis]) {\
  132.   double temp = min_array[axis]; min_array[axis] = max_array[axis]; max_array[axis] = temp;\
  133.   reverse_range[axis] = 1; \
  134.  } else reverse_range[axis] = (range_flags[axis]&RANGE_REVERSE); \
  135. }while(0)
  136.  
  137.  
  138. /* get optional [min:max] */
  139. #define LOAD_RANGE(axis) \
  140. do {\
  141.  if (equals(c_token, "[")) { \
  142.   c_token++; \
  143.   auto_array[axis] = load_range(axis,&min_array[axis], &max_array[axis], auto_array[axis]);\
  144.   if (!equals(c_token, "]"))\
  145.    int_error("']' expected", c_token);\
  146.   c_token++;\
  147.  }\
  148. } while (0)
  149.  
  150.  
  151. /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
  152.  * Do OUT_ACTION or UNDEF_ACTION as appropriate
  153.  * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
  154.  * VALUE must not be same as STORE
  155.  */
  156.  
  157. #define STORE_WITH_LOG_AND_FIXUP_RANGE(STORE, VALUE, TYPE, AXIS, OUT_ACTION, UNDEF_ACTION)\
  158. do { if (log_array[AXIS]) { if (VALUE<0.0) {TYPE = UNDEFINED; UNDEF_ACTION; break;} \
  159.               else if (VALUE == 0.0){STORE = -VERYLARGE; TYPE = OUTRANGE; OUT_ACTION; break;} \
  160.               else { STORE = log(VALUE)/log_base_array[AXIS]; } \
  161.      } else STORE = VALUE; \
  162.      if (TYPE != INRANGE) break;  /* dont set y range if x is outrange, for example */ \
  163.      if ( VALUE<min_array[AXIS] ) { \
  164.       if (auto_array[AXIS] & 1) min_array[AXIS] = VALUE; else { TYPE = OUTRANGE; OUT_ACTION; break; }  \
  165.      } \
  166.      if ( VALUE>max_array[AXIS] ) { \
  167.       if (auto_array[AXIS] & 2) max_array[AXIS] = VALUE; else { TYPE = OUTRANGE; OUT_ACTION; }   \
  168.      } \
  169. } while(0)
  170.  
  171. /* use this instead empty macro arguments to work around NeXT cpp bug */
  172. /* if this fails on any system, we might use ((void)0) */
  173. #define NOOP            /* */
  174.  
  175. /* check range and take logs of min and max if logscale
  176.  * this also restores min and max for ranges like [10:-10]
  177.  */
  178.  
  179. #ifdef HAVE_STRINGIZE
  180. # define RANGE_MSG(x) #x " range is less than threshold : see `set zero`"
  181. #else
  182. # define RANGE_MSG(x) "x range is less than threshold : see `set zero`"
  183. #endif
  184.  
  185. #define FIXUP_RANGE_FOR_LOG(AXIS, WHICH) \
  186. do { if (reverse_range[AXIS]) { \
  187.       double temp = min_array[AXIS]; \
  188.       min_array[AXIS] = max_array[AXIS]; \
  189.       max_array[AXIS] = temp; \
  190.      }\
  191.      if (log_array[AXIS]) { \
  192.       if (min_array[AXIS] <= 0.0 || max_array[AXIS] <= 0.0) \
  193.        int_error(RANGE_MSG(WHICH), NO_CARET); \
  194.       min_array[AXIS] = log(min_array[AXIS])/log_base_array[AXIS]; \
  195.       max_array[AXIS] = log(max_array[AXIS])/log_base_array[AXIS];  \
  196.     } \
  197. } while(0)
  198.  
  199.  
  200.  
  201. /* support for dynamic size of input line */
  202.  
  203. void plot3drequest()
  204. /*
  205.  * in the parametric case we would say splot [u= -Pi:Pi] [v= 0:2*Pi] [-1:1]
  206.  * [-1:1] [-1:1] sin(v)*cos(u),sin(v)*cos(u),sin(u) in the non-parametric
  207.  * case we would say only splot [x= -2:2] [y= -5:5] sin(x)*cos(y)
  208.  * 
  209.  */
  210. {
  211.     TBOOLEAN changed;
  212.     int dummy_token0 = -1, dummy_token1 = -1;
  213.  
  214.     is_3d_plot = TRUE;
  215.  
  216.     if (parametric && strcmp(dummy_var[0], "t") == 0) {
  217.     strcpy(dummy_var[0], "u");
  218.     strcpy(dummy_var[1], "v");
  219.     }
  220.     autoscale_lx = autoscale_x;
  221.     autoscale_ly = autoscale_y;
  222.     autoscale_lz = autoscale_z;
  223.  
  224.     if (!term)            /* unknown */
  225.     int_error("use 'set term' to set terminal type first", c_token);
  226.  
  227.     if (equals(c_token, "[")) {
  228.     c_token++;
  229.     if (isletter(c_token)) {
  230.         if (equals(c_token + 1, "=")) {
  231.         dummy_token0 = c_token;
  232.         c_token += 2;
  233.         } else {
  234.         /* oops; probably an expression with a variable. */
  235.         /* Parse it as an xmin expression. */
  236.         /* used to be: int_error("'=' expected",c_token); */
  237.         }
  238.     }
  239.     changed = parametric ? load_range(U_AXIS, &umin, &umax, autoscale_lu) : load_range(FIRST_X_AXIS, &xmin, &xmax, autoscale_lx);
  240.     if (!equals(c_token, "]"))
  241.         int_error("']' expected", c_token);
  242.     c_token++;
  243.     /* if (changed) */
  244.     if (parametric)
  245.         /* autoscale_lu = FALSE; */
  246.         autoscale_lu = changed;
  247.     else
  248.         /* autoscale_lx = FALSE; */
  249.         autoscale_lx = changed;
  250.     }
  251.     if (equals(c_token, "[")) {
  252.     c_token++;
  253.     if (isletter(c_token)) {
  254.         if (equals(c_token + 1, "=")) {
  255.         dummy_token1 = c_token;
  256.         c_token += 2;
  257.         } else {
  258.         /* oops; probably an expression with a variable. */
  259.         /* Parse it as an xmin expression. */
  260.         /* used to be: int_error("'=' expected",c_token); */
  261.         }
  262.     }
  263.     changed = parametric ? load_range(V_AXIS, &vmin, &vmax, autoscale_lv) : load_range(FIRST_Y_AXIS, &ymin, &ymax, autoscale_ly);
  264.     if (!equals(c_token, "]"))
  265.         int_error("']' expected", c_token);
  266.     c_token++;
  267.     /* if (changed) */
  268.     if (parametric)
  269.         /* autoscale_lv = FALSE; */
  270.         autoscale_lv = changed;
  271.     else
  272.         /* autoscale_ly = FALSE; */
  273.         autoscale_ly = changed;
  274.     }
  275.     if (parametric) {
  276.     /* set optional x (parametric) or z ranges */
  277.     if (equals(c_token, "[")) {
  278.         c_token++;
  279.         autoscale_lx = load_range(FIRST_X_AXIS, &xmin, &xmax, autoscale_lx);
  280.         if (!equals(c_token, "]"))
  281.         int_error("']' expected", c_token);
  282.         c_token++;
  283.     }
  284.     /* set optional y ranges */
  285.     if (equals(c_token, "[")) {
  286.         c_token++;
  287.         autoscale_ly = load_range(FIRST_Y_AXIS, &ymin, &ymax, autoscale_ly);
  288.         if (!equals(c_token, "]"))
  289.         int_error("']' expected", c_token);
  290.         c_token++;
  291.     }
  292.     }                /* parametric */
  293.     if (equals(c_token, "[")) {    /* set optional z ranges */
  294.     c_token++;
  295.     autoscale_lz = load_range(FIRST_Z_AXIS, &zmin, &zmax, autoscale_lz);
  296.     if (!equals(c_token, "]"))
  297.         int_error("']' expected", c_token);
  298.     c_token++;
  299.     }
  300.     CHECK_REVERSE(FIRST_X_AXIS);
  301.     CHECK_REVERSE(FIRST_Y_AXIS);
  302.     CHECK_REVERSE(FIRST_Z_AXIS);
  303.  
  304.     /* use the default dummy variable unless changed */
  305.     if (dummy_token0 >= 0)
  306.     copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);
  307.     else
  308.     (void) strcpy(c_dummy_var[0], dummy_var[0]);
  309.  
  310.     if (dummy_token1 >= 0)
  311.     copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);
  312.     else
  313.     (void) strcpy(c_dummy_var[1], dummy_var[1]);
  314.  
  315.     eval_3dplots();
  316. }
  317.  
  318.  
  319. static void grid_nongrid_data(this_plot)
  320. struct surface_points *this_plot;
  321. {
  322.     int i, j, k;
  323.     double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;
  324.     struct iso_curve *old_iso_crvs = this_plot->iso_crvs;
  325.     struct iso_curve *icrv, *oicrv, *oicrvs;
  326.  
  327.     /* Compute XY bounding box on the original data. */
  328.     xmin = xmax = old_iso_crvs->points[0].x;
  329.     ymin = ymax = old_iso_crvs->points[0].y;
  330.     for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {
  331.     struct coordinate GPHUGE *points = icrv->points;
  332.  
  333.     for (i = 0; i < icrv->p_count; i++, points++) {
  334.         if (xmin > points->x)
  335.         xmin = points->x;
  336.         if (xmax < points->x)
  337.         xmax = points->x;
  338.         if (ymin > points->y)
  339.         ymin = points->y;
  340.         if (ymax < points->y)
  341.         ymax = points->y;
  342.     }
  343.     }
  344.  
  345.     dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);
  346.     dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);
  347.  
  348.     /* Create the new grid structure, and compute the low pass filtering from
  349.      * non grid to grid structure.
  350.      */
  351.     this_plot->iso_crvs = NULL;
  352.     this_plot->num_iso_read = dgrid3d_col_fineness;
  353.     this_plot->has_grid_topology = TRUE;
  354.     for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {
  355.     struct coordinate GPHUGE *points;
  356.  
  357.     icrv = iso_alloc(dgrid3d_row_fineness + 1);
  358.     icrv->p_count = dgrid3d_row_fineness;
  359.     icrv->next = this_plot->iso_crvs;
  360.     this_plot->iso_crvs = icrv;
  361.     points = icrv->points;
  362.  
  363.     for (j = 0, y = ymin; j < dgrid3d_row_fineness; j++, y += dy, points++) {
  364.         z = w = 0.0;
  365.  
  366. #ifndef BUGGY_DGRID_RANGING /* HBB 981209 */
  367.         /* as soon as ->type is changed to UNDEFINED, break out of
  368.          * two inner loops! */
  369.         points->type = INRANGE;
  370. #endif
  371.         for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
  372.         struct coordinate GPHUGE *opoints = oicrv->points;
  373.         for (k = 0; k < oicrv->p_count; k++, opoints++) {
  374.             double dist, dist_x = fabs(opoints->x - x), dist_y = fabs(opoints->y - y);
  375.  
  376.             switch (dgrid3d_norm_value) {
  377.             case 1:
  378.             dist = dist_x + dist_y;
  379.             break;
  380.             case 2:
  381.             dist = dist_x * dist_x + dist_y * dist_y;
  382.             break;
  383.             case 4:
  384.             dist = dist_x * dist_x + dist_y * dist_y;
  385.             dist *= dist;
  386.             break;
  387.             case 8:
  388.             dist = dist_x * dist_x + dist_y * dist_y;
  389.             dist *= dist;
  390.             dist *= dist;
  391.             break;
  392.             case 16:
  393.             dist = dist_x * dist_x + dist_y * dist_y;
  394.             dist *= dist;
  395.             dist *= dist;
  396.             dist *= dist;
  397.             break;
  398.             default:
  399.             dist = pow(dist_x, (double) dgrid3d_norm_value) +
  400.                 pow(dist_y, (double) dgrid3d_norm_value);
  401.             break;
  402.             }
  403.  
  404.             /* The weight of this point is inverse proportional
  405.              * to the distance.
  406.              */
  407.             if (dist == 0.0) {
  408. #ifndef BUGGY_DGRID_RANGING
  409.                 /* HBB 981209: revised flagging as undefined */
  410.             /* Supporting all those infinities on various
  411.              * platforms becomes tiresome, to say the least :-(
  412.              * Let's just return the first z where this happens,
  413.              * unchanged, and be done with this, period. */
  414.                 points->type = UNDEFINED;
  415.             z = opoints->z;
  416.             w = 1.0;
  417.             break; /* out of for (k...) loop */
  418. #else
  419. #if !defined(AMIGA_SC_6_1) && !defined(__PUREC__)
  420.             dist = VERYLARGE;
  421. #else /* !AMIGA_SC_6_1 && !__PUREC__ */
  422.             /* Multiplying VERYLARGE by opoints->z below
  423.              * might yield Inf (i.e. a number that can't
  424.              * be represented on the machine). This will
  425.              * result in points->z being set to NaN. It's
  426.              * better to have a pretty large number that is
  427.              * also on the safe side... The numbers that are
  428.              * read by gnuplot are float values anyway, so
  429.              * they can't be bigger than FLT_MAX. So setting
  430.              * dist to FLT_MAX^2 will make dist pretty large
  431.              * with respect to any value that has been read. */
  432.             dist = ((double) FLT_MAX) * ((double) FLT_MAX);
  433. #endif /* !AMIGA_SC_6_1 && !__PUREC__ */
  434. #endif /* BUGGY_DGRID_RANGING */
  435.             } else
  436.             dist = 1.0 / dist;
  437.  
  438.             z += opoints->z * dist;
  439.             w += dist;
  440.         }
  441. #ifndef BUGGY_DGRID_RANGING
  442.         if (points->type != INRANGE)
  443.           break; /* out of the second-inner loop as well ... */
  444. #endif
  445.         }
  446.  
  447. #ifndef BUGGY_DGRID_RANGING
  448.         /* Now that we've escaped the loops safely, we know that we
  449.          * do have a good value in z and w, so we can proceed just as
  450.          * if nothing had happened at all. Nice, isn't it? */
  451.         points->type = INRANGE;
  452.  
  453.         STORE_WITH_LOG_AND_FIXUP_RANGE(points->x, x, points->type, x_axis, NOOP, continue);
  454.         STORE_WITH_LOG_AND_FIXUP_RANGE(points->y, y, points->type, y_axis, NOOP, continue);
  455.         STORE_WITH_LOG_AND_FIXUP_RANGE(points->z, z / w, points->type, z_axis, NOOP, continue);
  456. #else
  457.         /* HBB 981026: original, short version of this code */
  458.         points->x = x;
  459.         points->y = y;
  460.         points->z = z / w;
  461.         points->type = INRANGE;
  462. #endif
  463.     }
  464.     }
  465.  
  466.     /* Delete the old non grid data. */
  467.     for (oicrvs = old_iso_crvs; oicrvs != NULL;) {
  468.     oicrv = oicrvs;
  469.     oicrvs = oicrvs->next;
  470.     iso_free(oicrv);
  471.     }
  472. }
  473.  
  474. static void get_3ddata(this_plot)
  475. struct surface_points *this_plot;
  476. /* this_plot->token is end of datafile spec, before title etc
  477.  * will be moved passed title etc after we return
  478.  */
  479. {
  480.     int xdatum = 0;
  481.     int ydatum = 0;
  482.     int i, j;
  483.     double v[3];
  484.     int pt_in_iso_crv = 0;
  485.     struct iso_curve *this_iso;
  486.  
  487.     if (mapping3d == MAP3D_CARTESIAN) {
  488.     if (df_no_use_specs == 2)
  489.         int_error("Need 1 or 3 columns for cartesian data", this_plot->token);
  490.     } else {
  491.     if (df_no_use_specs == 1)
  492.         int_error("Need 2 or 3 columns for polar data", this_plot->token);
  493.     }
  494.  
  495.     this_plot->num_iso_read = 0;
  496.     this_plot->has_grid_topology = TRUE;
  497.  
  498.     /* we ought to keep old memory - most likely case
  499.      * is a replot, so it will probably exactly fit into
  500.      * memory already allocated ?
  501.      */
  502.     if (this_plot->iso_crvs != NULL) {
  503.     struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;
  504.  
  505.     while (icrvs) {
  506.         icrv = icrvs;
  507.         icrvs = icrvs->next;
  508.         iso_free(icrv);
  509.     }
  510.     this_plot->iso_crvs = NULL;
  511.     }
  512.     /* data file is already open */
  513.  
  514.     if (df_matrix)
  515.     xdatum = df_3dmatrix(this_plot);
  516.     else {
  517.     /*{{{  read surface from text file */
  518.     struct iso_curve *this_iso = iso_alloc(samples);
  519.     struct coordinate GPHUGE *cp;
  520.     double x, y, z;
  521.  
  522.     while ((j = df_readline(v, 3)) != DF_EOF) {
  523.         if (j == DF_SECOND_BLANK)
  524.         break;        /* two blank lines */
  525.         if (j == DF_FIRST_BLANK) {
  526.         /* one blank line */
  527.         if (pt_in_iso_crv == 0) {
  528.             if (xdatum == 0)
  529.             continue;
  530.             pt_in_iso_crv = xdatum;
  531.         }
  532.         if (xdatum > 0) {
  533.             this_iso->p_count = xdatum;
  534.             this_iso->next = this_plot->iso_crvs;
  535.             this_plot->iso_crvs = this_iso;
  536.             this_plot->num_iso_read++;
  537.  
  538.             if (xdatum != pt_in_iso_crv)
  539.             this_plot->has_grid_topology = FALSE;
  540.  
  541.             this_iso = iso_alloc(pt_in_iso_crv);
  542.             xdatum = 0;
  543.             ydatum++;
  544.         }
  545.         continue;
  546.         }
  547.         /* its a data point or undefined */
  548.  
  549.         if (xdatum >= this_iso->p_max) {
  550.         /*
  551.          * overflow about to occur. Extend size of points[] array. We
  552.          * either double the size, or add 1000 points, whichever is a
  553.          * smaller increment. Note i = p_max.
  554.          */
  555.         iso_extend(this_iso,
  556.                xdatum + (xdatum < 1000 ? xdatum : 1000));
  557.         }
  558.         cp = this_iso->points + xdatum;
  559.  
  560.         if (j == DF_UNDEFINED) {
  561.         cp->type = UNDEFINED;
  562.         continue;
  563.         }
  564.         cp->type = INRANGE;    /* unless we find out different */
  565.  
  566.         switch (mapping3d) {
  567.         case MAP3D_CARTESIAN:
  568.         switch (j) {
  569.         case 1:
  570.             x = xdatum;
  571.             y = ydatum;
  572.             z = v[0];
  573.             break;
  574.         case 3:
  575.             x = v[0];
  576.             y = v[1];
  577.             z = v[2];
  578.             break;
  579.         default:
  580.             {
  581.             char msg[80];
  582.             sprintf(msg, "Need 1 or 3 columns - line %d", df_line_number);
  583.             int_error(msg, this_plot->token);
  584.             return;    /* avoid gcc -Wuninitialised for x,y,z */
  585.             }
  586.         }
  587.         break;
  588.         case MAP3D_SPHERICAL:
  589.         if (j < 2)
  590.             int_error("Need 2 or 3 columns", this_plot->token);
  591.         if (j < 3)
  592.             v[2] = 1;    /* default radius */
  593.         if (angles_format == ANGLES_DEGREES) {
  594.             v[0] *= DEG2RAD;    /* Convert to radians. */
  595.             v[1] *= DEG2RAD;
  596.         }
  597.         x = v[2] * cos(v[0]) * cos(v[1]);
  598.         y = v[2] * sin(v[0]) * cos(v[1]);
  599.         z = v[2] * sin(v[1]);
  600.         break;
  601.         case MAP3D_CYLINDRICAL:
  602.         if (j < 2)
  603.             int_error("Need 2 or 3 columns", this_plot->token);
  604.         if (j < 3)
  605.             v[2] = 1;    /* default radius */
  606.         if (angles_format == ANGLES_DEGREES) {
  607.             v[0] *= DEG2RAD;    /* Convert to radians. */
  608.         }
  609.         x = v[2] * cos(v[0]);
  610.         y = v[2] * sin(v[0]);
  611.         z = v[1];
  612.         break;
  613.         default:
  614.         int_error("Internal error : Unknown mapping type", NO_CARET);
  615.         return;
  616.         }
  617.  
  618.         /* adjust for logscales. Set min/max and point types.
  619.          * store in cp
  620.          */
  621.         cp->type = INRANGE;
  622.         /* cannot use continue, as macro is wrapped in a loop.
  623.          * I regard this as correct goto use
  624.          */
  625.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->x, x, cp->type, x_axis, NOOP, goto come_here_if_undefined);
  626.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->y, y, cp->type, y_axis, NOOP, goto come_here_if_undefined);
  627.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->z, z, cp->type, z_axis, NOOP, goto come_here_if_undefined);
  628.  
  629.         /* some may complain, but I regard this as the correct use
  630.          * of goto
  631.          */
  632.       come_here_if_undefined:
  633.         ++xdatum;
  634.     }            /* end of whileloop - end of surface */
  635.  
  636.     if (xdatum > 0) {
  637.         this_plot->num_iso_read++;    /* Update last iso. */
  638.         this_iso->p_count = xdatum;
  639.  
  640.         this_iso->next = this_plot->iso_crvs;
  641.         this_plot->iso_crvs = this_iso;
  642.  
  643.         if (xdatum != pt_in_iso_crv)
  644.         this_plot->has_grid_topology = FALSE;
  645.  
  646.     } else {
  647.         iso_free(this_iso);    /* Free last allocation. */
  648.     }
  649.  
  650.     if (dgrid3d && this_plot->num_iso_read > 0)
  651.         grid_nongrid_data(this_plot);
  652.     /*}}} */
  653.     }
  654.  
  655.     if (this_plot->num_iso_read <= 1)
  656.     this_plot->has_grid_topology = FALSE;
  657.     if (this_plot->has_grid_topology && !hidden3d) {
  658.     struct iso_curve *new_icrvs = NULL;
  659.     int num_new_iso = this_plot->iso_crvs->p_count, len_new_iso = this_plot->num_iso_read;
  660.  
  661.     /* Now we need to set the other direction (pseudo) isolines. */
  662.     for (i = 0; i < num_new_iso; i++) {
  663.         struct iso_curve *new_icrv = iso_alloc(len_new_iso);
  664.  
  665.         new_icrv->p_count = len_new_iso;
  666.  
  667.         for (j = 0, this_iso = this_plot->iso_crvs;
  668.          this_iso != NULL;
  669.          j++, this_iso = this_iso->next) {
  670.         /* copy whole point struct to get type too.
  671.          * wasteful for windows, with padding */
  672.         /* more efficient would be extra pointer to same struct */
  673.         new_icrv->points[j] = this_iso->points[i];
  674.         }
  675.  
  676.         new_icrv->next = new_icrvs;
  677.         new_icrvs = new_icrv;
  678.     }
  679.  
  680.     /* Append the new iso curves after the read ones. */
  681.     for (this_iso = this_plot->iso_crvs;
  682.          this_iso->next != NULL;
  683.          this_iso = this_iso->next);
  684.     this_iso->next = new_icrvs;
  685.     }
  686. }
  687.  
  688.  
  689.  
  690. static void print_3dtable(pcount)
  691. int pcount;
  692. {
  693.     register struct surface_points *this_plot;
  694.     int i, curve, surface;
  695.     struct iso_curve *icrvs;
  696.     struct coordinate GPHUGE *points;
  697.  
  698.     for (surface = 0, this_plot = first_3dplot; surface < pcount;
  699.      this_plot = this_plot->next_sp, surface++) {
  700.     fprintf(gpoutfile, "\n#Surface %d of %d surfaces\n", surface, pcount);
  701.     icrvs = this_plot->iso_crvs;
  702.     curve = 0;
  703.  
  704.     if (draw_surface) {
  705.         /* only the curves in one direction */
  706.         while (icrvs && curve < this_plot->num_iso_read) {
  707.         fprintf(gpoutfile, "\n#IsoCurve %d, %d points\n#x y z type\n",
  708.             curve, icrvs->p_count);
  709.         for (i = 0, points = icrvs->points; i < icrvs->p_count; i++) {
  710.             fprintf(gpoutfile, "%g %g %g %c\n",
  711.                 points[i].x,
  712.                 points[i].y,
  713.                 points[i].z,
  714.                 points[i].type == INRANGE ? 'i'
  715.                 : points[i].type == OUTRANGE ? 'o'
  716.                 : 'u');
  717.         }
  718.         icrvs = icrvs->next;
  719.         curve++;
  720.         }
  721.         putc('\n', gpoutfile);
  722.     }
  723.     if (draw_contour) {
  724.         int number = 0;
  725.         struct gnuplot_contours *c = this_plot->contours;
  726.         while (c) {
  727.         int count = c->num_pts;
  728.         struct coordinate GPHUGE *p = c->coords;
  729.         if (c->isNewLevel)
  730.             /* dont display count - contour split across chunks */
  731.             /* put # in case user wants to use it for a plot */
  732.             /* double blank line to allow plot ... index ... */
  733.             fprintf(gpoutfile, "\n# Contour %d, label: %s\n", number++, c->label);
  734.         for (; --count >= 0; ++p)
  735.             fprintf(gpoutfile, "%g %g %g\n", p->x, p->y, p->z);
  736.         /* blank line between segments of same contour */
  737.         putc('\n', gpoutfile);
  738.         c = c->next;
  739.         }
  740.     }
  741.     }
  742.     fflush(gpoutfile);
  743. }
  744.  
  745.  
  746.  
  747. #define SET_DUMMY_RANGE(AXIS) \
  748. do{\
  749.  if (parametric || polar) { \
  750.   t_min = tmin; t_max = tmax;\
  751.  } else if (log_array[AXIS]) {\
  752.   if (min_array[AXIS] <= 0.0 || max_array[AXIS] <= 0.0)\
  753.    int_error("x/x2 range must be greater than 0 for log scale!", NO_CARET);\
  754.   t_min = log(min_array[AXIS])/log_base_array[AXIS]; t_max = log(max_array[AXIS])/log_base_array[AXIS];\
  755.  } else {\
  756.   t_min = min_array[AXIS]; t_max = max_array[AXIS];\
  757.  }\
  758.  t_step = (t_max - t_min) / (samples - 1); \
  759. }while(0)
  760.  
  761.  
  762. /*
  763.  * This parses the splot command after any range specifications. To support
  764.  * autoscaling on the x/z axis, we want any data files to define the x/y
  765.  * range, then to plot any functions using that range. We thus parse the
  766.  * input twice, once to pick up the data files, and again to pick up the
  767.  * functions. Definitions are processed twice, but that won't hurt.
  768.  * div - okay, it doesn't hurt, but every time an option as added for
  769.  * datafiles, code to parse it has to be added here. Change so that
  770.  * we store starting-token in the plot structure.
  771.  */
  772. static void eval_3dplots()
  773. {
  774.     int i, j;
  775.     struct surface_points **tp_3d_ptr;
  776.     int start_token, end_token;
  777.     int begin_token;
  778.     TBOOLEAN some_data_files = FALSE, some_functions = FALSE;
  779.     int plot_num, line_num, point_num, crnt_param = 0;    /* 0 = z, 1 = y, 2 = x */
  780.     char *xtitle;
  781.     char *ytitle;
  782.  
  783.     /* Reset first_3dplot. This is usually done at the end of this function.
  784.      * If there is an error within this function, the memory is left allocated,
  785.      * since we cannot call sp_free if the list is incomplete
  786.      */
  787.     first_3dplot = NULL;
  788.  
  789.     /* put stuff into arrays to simplify access */
  790.     INIT_ARRAYS(FIRST_X_AXIS, xmin, xmax, autoscale_lx, is_log_x, base_log_x, log_base_log_x, 0);
  791.     INIT_ARRAYS(FIRST_Y_AXIS, ymin, ymax, autoscale_ly, is_log_y, base_log_y, log_base_log_y, 0);
  792.     INIT_ARRAYS(FIRST_Z_AXIS, zmin, zmax, autoscale_lz, is_log_z, base_log_z, log_base_log_z, 1);
  793.  
  794.     x_axis = FIRST_X_AXIS;
  795.     y_axis = FIRST_Y_AXIS;
  796.     z_axis = FIRST_Z_AXIS;
  797.  
  798.     tp_3d_ptr = &(first_3dplot);
  799.     plot_num = 0;
  800.     line_num = 0;        /* default line type */
  801.     point_num = 0;        /* default point type */
  802.  
  803.     xtitle = NULL;
  804.     ytitle = NULL;
  805.  
  806.     begin_token = c_token;
  807.  
  808. /*** First Pass: Read through data files ***/
  809.     /*
  810.      * This pass serves to set the x/yranges and to parse the command, as
  811.      * well as filling in every thing except the function data. That is done
  812.      * after the x/yrange is defined.
  813.      */
  814.     while (TRUE) {
  815.     if (END_OF_COMMAND)
  816.         int_error("function to plt3d expected", c_token);
  817.  
  818.     start_token = c_token;
  819.  
  820.     if (is_definition(c_token)) {
  821.         define();
  822.     } else {
  823.         int specs;
  824.         struct surface_points *this_plot;
  825.  
  826.         if (isstring(c_token)) {    /* data file to plot */
  827.  
  828.         /*{{{  data file */
  829.         if (parametric && crnt_param != 0)
  830.             int_error("previous parametric function not fully specified", c_token);
  831.  
  832.         if (!some_data_files) {
  833.             if (autoscale_lx & 1) {
  834.             min_array[FIRST_X_AXIS] = VERYLARGE;
  835.             }
  836.             if (autoscale_lx & 2) {
  837.             max_array[FIRST_X_AXIS] = -VERYLARGE;
  838.             }
  839.             if (autoscale_ly & 1) {
  840.             min_array[FIRST_Y_AXIS] = VERYLARGE;
  841.             }
  842.             if (autoscale_ly & 2) {
  843.             max_array[FIRST_Y_AXIS] = -VERYLARGE;
  844.             }
  845.             some_data_files = TRUE;
  846.         }
  847.         if (*tp_3d_ptr)
  848.             this_plot = *tp_3d_ptr;
  849.         else {        /* no memory malloc()'d there yet */
  850.             /* Allocate enough isosamples and samples */
  851.             this_plot = sp_alloc(0, 0, 0, 0);
  852.             *tp_3d_ptr = this_plot;
  853.         }
  854.  
  855.         this_plot->plot_type = DATA3D;
  856.         this_plot->plot_style = data_style;
  857.  
  858.         specs = df_open(3);
  859.         /* parses all datafile-specific modifiers */
  860.         /* we will load the data after parsing title,with,... */
  861.  
  862.         /* for capture to key */
  863.         this_plot->token = end_token = c_token - 1;
  864.  
  865.         /* this_plot->token is temporary, for errors in get_3ddata() */
  866.  
  867.         if (datatype[FIRST_X_AXIS] == TIME) {
  868.             if (specs < 3)
  869.             int_error("Need full using spec for x time data",
  870.                   c_token);
  871.             df_timecol[0] = 1;
  872.         }
  873.         if (datatype[FIRST_Y_AXIS] == TIME) {
  874.             if (specs < 3)
  875.             int_error("Need full using spec for y time data",
  876.                   c_token);
  877.             df_timecol[1] = 1;
  878.         }
  879.         if (datatype[FIRST_Z_AXIS] == TIME) {
  880.             if (specs < 3)
  881.             df_timecol[0] = 1;
  882.             else
  883.             df_timecol[2] = 1;
  884.         }
  885.         /*}}} */
  886.  
  887.         } else {        /* function to plot */
  888.  
  889.         /*{{{  function */
  890.         ++plot_num;
  891.         if (parametric)    {
  892.             /* Rotate between x/y/z axes */
  893.             /* +2 same as -1, but beats -ve problem */
  894.             crnt_param = (crnt_param + 2) % 3;
  895.         }
  896.  
  897.         if (*tp_3d_ptr) {
  898.             this_plot = *tp_3d_ptr;
  899.             if (!hidden3d)
  900.             sp_replace(this_plot, samples_1, iso_samples_1,
  901.                    samples_2, iso_samples_2);
  902.             else
  903.             sp_replace(this_plot, iso_samples_1, 0,
  904.                    0, iso_samples_2);
  905.  
  906.         } else {    /* no memory malloc()'d there yet */
  907.             /* Allocate enough isosamples and samples */
  908.             if (!hidden3d)
  909.             this_plot = sp_alloc(samples_1, iso_samples_1,
  910.                          samples_2, iso_samples_2);
  911.             else
  912.             this_plot = sp_alloc(iso_samples_1, 0,
  913.                          0, iso_samples_2);
  914.             *tp_3d_ptr = this_plot;
  915.         }
  916.  
  917.         this_plot->plot_type = FUNC3D;
  918.         this_plot->has_grid_topology = TRUE;
  919.         this_plot->plot_style = func_style;
  920.         this_plot->num_iso_read = iso_samples_2;
  921.         dummy_func = &plot_func;
  922.         plot_func.at = temp_at();
  923.         dummy_func = NULL;
  924.         /* ignore it for now */
  925.         some_functions = TRUE;
  926.         end_token = c_token - 1;
  927.         /*}}} */
  928.  
  929.         }            /* end of IS THIS A FILE OR A FUNC block */
  930.  
  931.  
  932.         /*{{{  title */
  933.         if (this_plot->title) {
  934.         free(this_plot->title);
  935.         this_plot->title = NULL;
  936.         }
  937.         if (almost_equals(c_token, "t$itle")) {
  938.         /*                      if (!isstring(++c_token))
  939.            int_error("Expected title", c_token);
  940.            m_quote_capture(&(this_plot->title), c_token, c_token);
  941.          */
  942.         if (parametric) {
  943.             if (crnt_param != 0)
  944.             int_error("\"title\" allowed only after parametric function fully specified", c_token);
  945.             else {
  946.             if (xtitle != NULL)
  947.                 xtitle[0] = NUL;    /* Remove default title . */
  948.             if (ytitle != NULL)
  949.                 ytitle[0] = NUL;    /* Remove default title . */
  950.             }
  951.         }
  952.         if (isstring(++c_token))
  953.             m_quote_capture(&(this_plot->title), c_token, c_token);
  954.         else
  955.             int_error("expecting \"title\" for plot", c_token);
  956.         /* end of new method */
  957.         ++c_token;
  958.         } else if (almost_equals(c_token, "not$itle")) {
  959.         if (xtitle != NULL)
  960.             xtitle[0] = '\0';
  961.         if (ytitle != NULL)
  962.             ytitle[0] = '\0';
  963.         /*   this_plot->title = NULL;   */
  964.         ++c_token;
  965.         } else {
  966.         m_capture(&(this_plot->title), start_token, end_token);
  967.         if (crnt_param == 2)
  968.             xtitle = this_plot->title;
  969.         else if (crnt_param == 1)
  970.             ytitle = this_plot->title;
  971.         }
  972.         /*}}} */
  973.  
  974.  
  975.         /*{{{  line types, widths, ... */
  976.         this_plot->lp_properties.l_type = line_num;
  977.         this_plot->lp_properties.p_type = point_num;
  978.  
  979.         if (almost_equals(c_token, "w$ith")) {
  980.         this_plot->plot_style = get_style();
  981.         }
  982.         /* pick up line/point specs
  983.          * - point spec allowed if style uses points, ie style&2 != 0
  984.          * - keywords are optional
  985.          */
  986.         LP_PARSE(this_plot->lp_properties, 1, this_plot->plot_style & 2,
  987.              line_num, point_num);
  988.  
  989.         /* allow old-style syntax too - ignore case lt 3 4 for example */
  990.         if (!equals(c_token, ",") && !END_OF_COMMAND) {
  991.         struct value t;
  992.         this_plot->lp_properties.l_type =
  993.             this_plot->lp_properties.p_type = (int) real(const_express(&t)) - 1;
  994.  
  995.         if (!equals(c_token, ",") && !END_OF_COMMAND)
  996.             this_plot->lp_properties.p_type = (int) real(const_express(&t)) - 1;
  997.         }
  998.         if (this_plot->plot_style & 2)    /* lines, linesp, ... */
  999.         if (crnt_param == 0)
  1000.             point_num +=
  1001.             1 + (draw_contour != 0)
  1002.             + (hidden3d != 0);
  1003.  
  1004.         if (crnt_param == 0)
  1005.         line_num += 1 + (draw_contour != 0)
  1006.             + (hidden3d != 0);
  1007.         /*}}} */
  1008.  
  1009.  
  1010.         /* now get the data... having to think hard here...
  1011.          * first time through, we fill in this_plot. For second
  1012.          * surface in file, we have to allocate another surface
  1013.          * struct. BUT we may allocate this store only to
  1014.          * find that it is merely some blank lines at end of file
  1015.          * tp_3d_ptr is still pointing at next field of prev. plot,
  1016.          * before :    prev_or_first -> this_plot -> possible_preallocated_store
  1017.          *                tp_3d_ptr--^
  1018.          * after  :    prev_or_first -> first -> second -> last -> possibly_more_store
  1019.          *                                        tp_3d_ptr ----^
  1020.          * if file is empty, tp_3d_ptr is not moved. this_plot continues
  1021.          * to point at allocated storage, but that will be reused later
  1022.          */
  1023.  
  1024.         assert(this_plot == *tp_3d_ptr);
  1025.  
  1026.         if (this_plot->plot_type == DATA3D) {
  1027.         /*{{{  read data */
  1028.         /* remember settings for second surface in file */
  1029.         struct lp_style_type *these_props = &(this_plot->lp_properties);
  1030.         enum PLOT_STYLE this_style = this_plot->plot_style;
  1031.  
  1032.         int this_token = this_plot->token;
  1033.         while (!df_eof) {
  1034.             this_plot = *tp_3d_ptr;
  1035.             assert(this_plot != NULL);
  1036.  
  1037.             /* dont move tp_3d_ptr until we are sure we
  1038.              * have read a surface
  1039.              */
  1040.             
  1041.             /* used by get_3ddata() */
  1042.             this_plot->token = this_token;
  1043.  
  1044.             get_3ddata(this_plot);
  1045.             /* for second pass */
  1046.             this_plot->token = c_token;
  1047.  
  1048.             if (this_plot->num_iso_read == 0)
  1049.             /* probably df_eof, in which case we
  1050.              * will leave loop. if not eof, then
  1051.              * how come we got no surface ? - retry
  1052.              * in neither case do we update tp_3d_ptr
  1053.              */
  1054.             continue;
  1055.  
  1056.             /* okay, we have read a surface */
  1057.             ++plot_num;
  1058.             tp_3d_ptr = &(this_plot->next_sp);
  1059.             if (df_eof)
  1060.             break;
  1061.  
  1062.             /* there might be another surface so allocate
  1063.              * and prepare another surface structure
  1064.              * This does no harm if in fact there are
  1065.              * no more surfaces to read
  1066.              */
  1067.  
  1068.             if ((this_plot = *tp_3d_ptr) != NULL) {
  1069.             if (this_plot->title) {
  1070.                 free(this_plot->title);
  1071.                 this_plot->title = NULL;
  1072.             }
  1073.             } else {
  1074.             /* Allocate enough isosamples and samples */
  1075.             this_plot = *tp_3d_ptr = sp_alloc(0, 0, 0, 0);
  1076.             }
  1077.  
  1078.             this_plot->plot_type = DATA3D;
  1079.             this_plot->plot_style = this_style;
  1080.             /* Struct copy */
  1081.             this_plot->lp_properties = *these_props;
  1082.         }
  1083.         df_close();
  1084.         /*}}} */
  1085.  
  1086.         } else {        /* not a data file */
  1087.         tp_3d_ptr = &(this_plot->next_sp);
  1088.         this_plot->token = c_token;    /* store for second pass */
  1089.         }
  1090.  
  1091.     }        /* !is_definition() : end of scope of this_plot */
  1092.  
  1093.  
  1094.     if (equals(c_token, ","))
  1095.         c_token++;
  1096.     else
  1097.         break;
  1098.  
  1099.     }                /* while(TRUE), ie first pass */
  1100.  
  1101.  
  1102.     if (parametric && crnt_param != 0)
  1103.     int_error("parametric function not fully specified", NO_CARET);
  1104.  
  1105.  
  1106. /*** Second Pass: Evaluate the functions ***/
  1107.     /*
  1108.      * Everything is defined now, except the function data. We expect no
  1109.      * syntax errors, etc, since the above parsed it all. This makes the code
  1110.      * below simpler. If autoscale_ly, the yrange may still change.
  1111.      * - eh ?  - z can still change.  x/y/z can change if we are parametric ??
  1112.      */
  1113.  
  1114.     if (some_functions) {
  1115.  
  1116.     /* I've changed the controlled variable in fn plots to u_min etc since
  1117.      * it's easier for me to think parametric - 'normal' plot is after all
  1118.      * a special case. I was confused about x_min being both minimum of
  1119.      * x values found, and starting value for fn plots.
  1120.      */
  1121.     register double u_min, u_max, u_step, v_min, v_max, v_step;
  1122.     double uisodiff, visodiff;
  1123.     struct surface_points *this_plot;
  1124.  
  1125.  
  1126.     if (!parametric) {
  1127.         /*{{{  check ranges */
  1128.         /* give error if xrange badly set from missing datafile error
  1129.          * parametric fn can still set ranges
  1130.          * if there are no fns, we'll report it later as 'nothing to plot'
  1131.          */
  1132.  
  1133.         if (min_array[FIRST_X_AXIS] == VERYLARGE || 
  1134.         max_array[FIRST_X_AXIS] == -VERYLARGE) {
  1135.         int_error("x range is invalid", c_token);
  1136.         }
  1137.         if (min_array[FIRST_Y_AXIS] == VERYLARGE || 
  1138.         max_array[FIRST_Y_AXIS] == -VERYLARGE) {
  1139.         int_error("y range is invalid", c_token);
  1140.         }
  1141.         /* check that xmin -> xmax is not too small */
  1142.         fixup_range(FIRST_X_AXIS, "x");
  1143.         fixup_range(FIRST_Y_AXIS, "y");
  1144.         /*}}} */
  1145.     }
  1146.     if (parametric && !some_data_files) {
  1147.         /*{{{  set ranges */
  1148.         /* parametric fn can still change x/y range */
  1149.         if (autoscale_lx & 1)
  1150.         min_array[FIRST_X_AXIS] = VERYLARGE;
  1151.         if (autoscale_lx & 2)
  1152.         max_array[FIRST_X_AXIS] = -VERYLARGE;
  1153.         if (autoscale_ly & 1)
  1154.         min_array[FIRST_Y_AXIS] = VERYLARGE;
  1155.         if (autoscale_ly & 2)
  1156.         max_array[FIRST_Y_AXIS] = -VERYLARGE;
  1157.         /*}}} */
  1158.     }
  1159.     if (parametric) {
  1160.         u_min = umin;
  1161.         u_max = umax;
  1162.         v_min = vmin;
  1163.         v_max = vmax;
  1164.     } else {
  1165.         /*{{{  figure ranges, taking logs etc into account */
  1166.         if (is_log_x) {
  1167.         if (min_array[FIRST_X_AXIS] <= 0.0 ||
  1168.             max_array[FIRST_X_AXIS] <= 0.0)
  1169.             int_error("x range must be greater than 0 for log scale!",
  1170.                   NO_CARET);
  1171.         u_min = log(min_array[FIRST_X_AXIS]) / log_base_log_x;
  1172.         u_max = log(max_array[FIRST_X_AXIS]) / log_base_log_x;
  1173.         } else {
  1174.         u_min = min_array[FIRST_X_AXIS];
  1175.         u_max = max_array[FIRST_X_AXIS];
  1176.         }
  1177.  
  1178.         if (is_log_y) {
  1179.         if (min_array[FIRST_Y_AXIS] <= 0.0 ||
  1180.             max_array[FIRST_Y_AXIS] <= 0.0) {
  1181.             int_error("y range must be greater than 0 for log scale!",
  1182.                   NO_CARET);
  1183.         }
  1184.         v_min = log(min_array[FIRST_Y_AXIS]) / log_base_log_y;
  1185.         v_max = log(max_array[FIRST_Y_AXIS]) / log_base_log_y;
  1186.         } else {
  1187.         v_min = min_array[FIRST_Y_AXIS];
  1188.         v_max = max_array[FIRST_Y_AXIS];
  1189.         }
  1190.         /*}}} */
  1191.     }
  1192.  
  1193.  
  1194.     if (samples_1 < 2 || samples_2 < 2 || iso_samples_1 < 2 ||
  1195.         iso_samples_2 < 2) {
  1196.         int_error("samples or iso_samples < 2. Must be at least 2.",
  1197.               NO_CARET);
  1198.     }
  1199.  
  1200.     /* start over */
  1201.     this_plot = first_3dplot;
  1202.     c_token = begin_token;
  1203.  
  1204.     /* why do attributes of this_plot matter ? */
  1205.  
  1206.     if (this_plot && this_plot->has_grid_topology && hidden3d) {
  1207.         u_step = (u_max - u_min) / (iso_samples_1 - 1);
  1208.         v_step = (v_max - v_min) / (iso_samples_2 - 1);
  1209.     } else {
  1210.         u_step = (u_max - u_min) / (samples_1 - 1);
  1211.         v_step = (v_max - v_min) / (samples_2 - 1);
  1212.     }
  1213.  
  1214.     uisodiff = (u_max - u_min) / (iso_samples_1 - 1);
  1215.     visodiff = (v_max - v_min) / (iso_samples_2 - 1);
  1216.  
  1217.  
  1218.     /* Read through functions */
  1219.     while (TRUE) {
  1220.         if (is_definition(c_token)) {
  1221.         define();
  1222.         } else {
  1223.  
  1224.         if (!isstring(c_token)) {    /* func to plot */
  1225.             /*{{{  evaluate function */
  1226.             struct iso_curve *this_iso = this_plot->iso_crvs;
  1227.             struct coordinate GPHUGE *points = this_iso->points;
  1228.             int num_sam_to_use, num_iso_to_use;
  1229.  
  1230.             if (parametric)
  1231.             crnt_param = (crnt_param + 2) % 3;
  1232.  
  1233.             dummy_func = &plot_func;
  1234.             plot_func.at = temp_at();    /* reparse function */
  1235.             dummy_func = NULL;
  1236.             num_iso_to_use = iso_samples_2;
  1237.  
  1238.             if (!(this_plot->has_grid_topology && hidden3d))
  1239.             num_sam_to_use = samples_1;
  1240.             else
  1241.             num_sam_to_use = iso_samples_1;
  1242.  
  1243.             for (j = 0; j < num_iso_to_use; j++) {
  1244.             double y = v_min + j * visodiff;
  1245.             /* if (is_log_y) PEM fix logscale y axis */
  1246.             /* y = pow(log_base_log_y,y); 26-Sep-89 */
  1247.             /* parametric => NOT a log quantity (?) */
  1248.             (void) Gcomplex(&plot_func.dummy_values[1],
  1249.             !parametric && is_log_y ? pow(base_log_y, y) : y,
  1250.                     0.0);
  1251.  
  1252.             for (i = 0; i < num_sam_to_use; i++) {
  1253.                 double x = u_min + i * u_step;
  1254.                 struct value a;
  1255.                 double temp;
  1256.  
  1257.                 /* if (is_log_x) PEM fix logscale x axis */
  1258.                 /* x = pow(base_log_x,x); 26-Sep-89 */
  1259.                 /* parametric => NOT a log quantity (?) */
  1260.                 (void) Gcomplex(&plot_func.dummy_values[0],
  1261.                         !parametric && is_log_x ? pow(base_log_x, x) : x, 0.0);
  1262.  
  1263.                 points[i].x = x;
  1264.                 points[i].y = y;
  1265.  
  1266.                 evaluate_at(plot_func.at, &a);
  1267.  
  1268.                 if (undefined || (fabs(imag(&a)) > zero)) {
  1269.                 points[i].type = UNDEFINED;
  1270.                 continue;
  1271.                 }
  1272.                 temp = real(&a);
  1273.  
  1274.                 points[i].type = INRANGE;
  1275.                 STORE_WITH_LOG_AND_FIXUP_RANGE(points[i].z, temp, points[i].type, crnt_param, NOOP, NOOP);
  1276.  
  1277.             }
  1278.             this_iso->p_count = num_sam_to_use;
  1279.             this_iso = this_iso->next;
  1280.             points = this_iso ? this_iso->points : NULL;
  1281.             }
  1282.  
  1283.             if (!(this_plot->has_grid_topology && hidden3d)) {
  1284.             num_iso_to_use = iso_samples_1;
  1285.             num_sam_to_use = samples_2;
  1286.             for (i = 0; i < num_iso_to_use; i++) {
  1287.                 double x = u_min + i * uisodiff;
  1288.                 /* if (is_log_x) PEM fix logscale x axis */
  1289.                 /* x = pow(base_log_x,x); 26-Sep-89 */
  1290.                 /* if parametric, no logs involved - 3.6 */
  1291.                 (void) Gcomplex(&plot_func.dummy_values[0],
  1292.                         (!parametric && is_log_x) ? pow(base_log_x, x) : x, 0.0);
  1293.  
  1294.                 for (j = 0; j < num_sam_to_use; j++) {
  1295.                 double y = v_min + j * v_step;
  1296.                 struct value a;
  1297.                 double temp;
  1298.                 /* if (is_log_y) PEM fix logscale y axis */
  1299.                 /* y = pow(base_log_y,y); 26-Sep-89 */
  1300.                 (void) Gcomplex(&plot_func.dummy_values[1],
  1301.                         (!parametric && is_log_y) ? pow(base_log_y, y) : y, 0.0);
  1302.  
  1303.                 points[j].x = x;
  1304.                 points[j].y = y;
  1305.  
  1306.                 evaluate_at(plot_func.at, &a);
  1307.  
  1308.                 if (undefined || (fabs(imag(&a)) > zero)) {
  1309.                     points[j].type = UNDEFINED;
  1310.                     continue;
  1311.                 }
  1312.                 temp = real(&a);
  1313.  
  1314.                 points[j].type = INRANGE;
  1315.                 STORE_WITH_LOG_AND_FIXUP_RANGE(points[j].z, temp, points[j].type, crnt_param, NOOP, NOOP);
  1316.                 }
  1317.                 this_iso->p_count = num_sam_to_use;
  1318.                 this_iso = this_iso->next;
  1319.                 points = this_iso ? this_iso->points : NULL;
  1320.             }
  1321.             }
  1322.             /*}}} */
  1323.         }        /* end of ITS A FUNCTION TO PLOT */
  1324.  
  1325.         /* we saved it from first pass */ 
  1326.         c_token = this_plot->token;
  1327.  
  1328.         /* one data file can make several plots */
  1329.         do
  1330.             this_plot = this_plot->next_sp;
  1331.         while (this_plot && this_plot->token == c_token);
  1332.  
  1333.         }            /* !is_definition */
  1334.  
  1335.         if (equals(c_token, ","))
  1336.         c_token++;
  1337.         else
  1338.         break;
  1339.  
  1340.     }            /* while(TRUE) */
  1341.  
  1342.  
  1343.     if (parametric) {
  1344.         /* Now actually fix the plot triplets to be single plots. */
  1345.         parametric_3dfixup(first_3dplot, &plot_num);
  1346.     }
  1347.     }                /* some functions */
  1348.     /* if first_3dplot is NULL, we have no functions or data at all.
  1349.        * This can happen, if you type "splot x=5", since x=5 is a
  1350.        * variable assignment
  1351.      */
  1352.     if (plot_num == 0 || first_3dplot == NULL) {
  1353.     int_error("no functions or data to plot", c_token);
  1354.     }
  1355.     if (min_array[FIRST_X_AXIS] == VERYLARGE ||
  1356.     max_array[FIRST_X_AXIS] == -VERYLARGE ||
  1357.     min_array[FIRST_Y_AXIS] == VERYLARGE ||
  1358.     max_array[FIRST_Y_AXIS] == -VERYLARGE ||
  1359.     min_array[FIRST_Z_AXIS] == VERYLARGE ||
  1360.     max_array[FIRST_Z_AXIS] == -VERYLARGE)
  1361.     int_error("All points undefined", NO_CARET);
  1362.  
  1363.     fixup_range(FIRST_X_AXIS, "x");
  1364.     fixup_range(FIRST_Y_AXIS, "y");
  1365.     fixup_range(FIRST_Z_AXIS, "z");
  1366.  
  1367.     FIXUP_RANGE_FOR_LOG(FIRST_X_AXIS, x);
  1368.     FIXUP_RANGE_FOR_LOG(FIRST_Y_AXIS, y);
  1369.     FIXUP_RANGE_FOR_LOG(FIRST_Z_AXIS, z);
  1370.  
  1371.     /* last parameter should take plot size into effect...
  1372.      * probably needs to be moved to graph3d.c
  1373.      * in the meantime, a value of 20 gives same behaviour
  1374.      * as 3.5 which will do for the moment
  1375.      */
  1376.  
  1377.     if (xtics)
  1378.     setup_tics(FIRST_X_AXIS, &xticdef, xformat, 20);
  1379.     if (ytics)
  1380.     setup_tics(FIRST_Y_AXIS, &yticdef, yformat, 20);
  1381.     if (ztics)
  1382.     setup_tics(FIRST_Z_AXIS, &zticdef, zformat, 20);
  1383.  
  1384. #define WRITEBACK(axis,min,max) \
  1385. if(range_flags[axis]&RANGE_WRITEBACK) \
  1386.   {if (auto_array[axis]&1) min = min_array[axis]; \
  1387.    if (auto_array[axis]&2) max = max_array[axis]; \
  1388.   }
  1389.  
  1390.     WRITEBACK(FIRST_X_AXIS, xmin, xmax);
  1391.     WRITEBACK(FIRST_Y_AXIS, ymin, ymax);
  1392.     WRITEBACK(FIRST_Z_AXIS, zmin, zmax);
  1393.  
  1394.     if (plot_num == 0 || first_3dplot == NULL) {
  1395.         int_error("no functions or data to plot", c_token);
  1396.     }
  1397.     /* Creates contours if contours are to be plotted as well. */
  1398.  
  1399.     if (draw_contour) {
  1400.     struct surface_points *this_plot;
  1401.     for (this_plot = first_3dplot, i = 0;
  1402.          i < plot_num;
  1403.          this_plot = this_plot->next_sp, i++) {
  1404.  
  1405.         if (this_plot->contours) {
  1406.         struct gnuplot_contours *cntrs = this_plot->contours;
  1407.  
  1408.         while (cntrs) {
  1409.             struct gnuplot_contours *cntr = cntrs;
  1410.             cntrs = cntrs->next;
  1411.             free(cntr->coords);
  1412.             free(cntr);
  1413.         }
  1414.         }
  1415.         /* Make sure this one can be contoured. */
  1416.         if (!this_plot->has_grid_topology) {
  1417.         this_plot->contours = NULL;
  1418.         fputs("Notice: cannot contour non grid data!\n", stderr);
  1419.         /* changed from int_error by recommendation of
  1420.          * rkc@xn.ll.mit.edu
  1421.          */
  1422.         } else if (this_plot->plot_type == DATA3D) {
  1423.         this_plot->contours = contour(
  1424.                          this_plot->num_iso_read,
  1425.                          this_plot->iso_crvs,
  1426.                          contour_levels, contour_pts,
  1427.                          contour_kind, contour_order,
  1428.                            levels_kind, levels_list);
  1429.         } else {
  1430.         this_plot->contours = contour(iso_samples_2,
  1431.                           this_plot->iso_crvs,
  1432.                           contour_levels, contour_pts,
  1433.                           contour_kind, contour_order,
  1434.                           levels_kind, levels_list);
  1435.         }
  1436.     }
  1437.     }                /* draw_contour */
  1438.     /* perform the plot */
  1439.     if (strcmp(term->name, "table") == 0)
  1440.     print_3dtable(plot_num);
  1441.     else {
  1442.     START_LEAK_CHECK();    /* assert no memory leaks here ! */
  1443.     do_3dplot(first_3dplot, plot_num);
  1444.     END_LEAK_CHECK();
  1445.     }
  1446.  
  1447.     /* if we get here, all went well, so record the line for replot */
  1448.     if (plot_token != -1) {
  1449.     /* note that m_capture also frees the old replot_line */
  1450.     m_capture(&replot_line, plot_token, c_token - 1);
  1451.     plot_token = -1;
  1452.     }
  1453.     sp_free(first_3dplot);
  1454.     first_3dplot = NULL;
  1455. }
  1456.  
  1457.  
  1458.  
  1459.  
  1460.  
  1461. static void parametric_3dfixup(start_plot, plot_num)
  1462. struct surface_points *start_plot;
  1463. int *plot_num;
  1464. /*
  1465.  * The hardest part of this routine is collapsing the FUNC plot types in the
  1466.  * list (which are gauranteed to occur in (x,y,z) triplets while preserving
  1467.  * the non-FUNC type plots intact.  This means we have to work our way
  1468.  * through various lists.  Examples (hand checked):
  1469.  * start_plot:F1->F2->F3->NULL ==> F3->NULL
  1470.  * start_plot:F1->F2->F3->F4->F5->F6->NULL ==> F3->F6->NULL
  1471.  * start_plot:F1->F2->F3->D1->D2->F4->F5->F6->D3->NULL ==>
  1472.  * F3->D1->D2->F6->D3->NULL
  1473.  */
  1474. {
  1475. /*
  1476.  * I initialized *free_list with NULL, because my compiler warns some lines
  1477.  * later that it might be uninited. The code however seems to not access that
  1478.  * line in that case, but if I'm right, my change is OK and if not, this is a
  1479.  * serious bug in the code.
  1480.  *
  1481.  * x and y ranges now fixed in eval_3dplots
  1482.  */
  1483.     struct surface_points *xp, *new_list, *free_list = NULL;
  1484.     struct surface_points **last_pointer = &new_list;
  1485.  
  1486.     int i, tlen, surface;
  1487.     char *new_title;
  1488.  
  1489.     /*
  1490.      * Ok, go through all the plots and move FUNC3D types together.  Note:
  1491.      * this originally was written to look for a NULL next pointer, but
  1492.      * gnuplot wants to be sticky in grabbing memory and the right number of
  1493.      * items in the plot list is controlled by the plot_num variable.
  1494.      * 
  1495.      * Since gnuplot wants to do this sticky business, a free_list of
  1496.      * surface_points is kept and then tagged onto the end of the plot list
  1497.      * as this seems more in the spirit of the original memory behavior than
  1498.      * simply freeing the memory.  I'm personally not convinced this sort of
  1499.      * concern is worth it since the time spent computing points seems to
  1500.      * dominate any garbage collecting that might be saved here...
  1501.      */
  1502.     new_list = xp = start_plot;
  1503.     for (surface = 0; surface < *plot_num; surface++) {
  1504.     if (xp->plot_type == FUNC3D) {
  1505.         struct surface_points *yp = xp->next_sp;
  1506.         struct surface_points *zp = yp->next_sp;
  1507.  
  1508.         /* Here's a FUNC3D parametric function defined as three parts.
  1509.          * Go through all the points and assign the x's and y's from xp and
  1510.          * yp to zp. min/max already done
  1511.          */
  1512.         struct iso_curve *xicrvs = xp->iso_crvs;
  1513.         struct iso_curve *yicrvs = yp->iso_crvs;
  1514.         struct iso_curve *zicrvs = zp->iso_crvs;
  1515.  
  1516.         (*plot_num) -= 2;
  1517.  
  1518.         assert(INRANGE < OUTRANGE && OUTRANGE < UNDEFINED);
  1519.  
  1520.         while (zicrvs) {
  1521.         struct coordinate GPHUGE *xpoints = xicrvs->points,
  1522.          GPHUGE * ypoints = yicrvs->points, GPHUGE * zpoints = zicrvs->points;
  1523.         for (i = 0; i < zicrvs->p_count; ++i) {
  1524.             zpoints[i].x = xpoints[i].z;
  1525.             zpoints[i].y = ypoints[i].z;
  1526.             if (zpoints[i].type < xpoints[i].type)
  1527.             zpoints[i].type = xpoints[i].type;
  1528.             if (zpoints[i].type < ypoints[i].type)
  1529.             zpoints[i].type = ypoints[i].type;
  1530.  
  1531.         }
  1532.         xicrvs = xicrvs->next;
  1533.         yicrvs = yicrvs->next;
  1534.         zicrvs = zicrvs->next;
  1535.         }
  1536.  
  1537.         /* Ok, fix up the title to include xp and yp plots. */
  1538.         if (((xp->title && xp->title[0] != '\0') ||
  1539.          (yp->title && yp->title[0] != '\0')) && zp->title) {
  1540.         tlen = (xp->title ? strlen(xp->title) : 0) +
  1541.             (yp->title ? strlen(yp->title) : 0) +
  1542.             (zp->title ? strlen(zp->title) : 0) + 5;
  1543.         new_title = gp_alloc((unsigned long) tlen, "string");
  1544.         new_title[0] = 0;
  1545.         if (xp->title && xp->title[0] != '\0') {
  1546.             strcat(new_title, xp->title);
  1547.             strcat(new_title, ", ");    /* + 2 */
  1548.         }
  1549.         if (yp->title && yp->title[0] != '\0') {
  1550.             strcat(new_title, yp->title);
  1551.             strcat(new_title, ", ");    /* + 2 */
  1552.         }
  1553.         strcat(new_title, zp->title);
  1554.         free(zp->title);
  1555.         zp->title = new_title;
  1556.         }
  1557.         /* add xp and yp to head of free list */
  1558.         assert(xp->next_sp == yp);
  1559.         yp->next_sp = free_list;
  1560.         free_list = xp;
  1561.  
  1562.         /* add zp to tail of new_list */
  1563.         *last_pointer = zp;
  1564.         last_pointer = &(zp->next_sp);
  1565.         xp = zp->next_sp;
  1566.     } else {        /* its a data plot */
  1567.         assert(*last_pointer == xp);    /* think this is true ! */
  1568.         last_pointer = &(xp->next_sp);
  1569.         xp = xp->next_sp;
  1570.     }
  1571.     }
  1572.  
  1573.     /* Ok, append free list and write first_plot */
  1574.     *last_pointer = free_list;
  1575.     first_3dplot = new_list;
  1576. }
  1577.